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Abstract. Many physical phenomena, governed by partial differential equa- 
tions (PDEs), are second order in nature. This makes sense to pose the control 
on the second order derivatives of the field solution, in addition to zero and 

U first order ones, to consistently control the underlaying process. However, this 

type of control is nontrivial and to the best of our knowledge there is nigher a 
C 3 ■ theoretic nor a numeric work in this regard. The present work goals to do the 

first quest in this regard, examining a problem of this type using a numerical 
simulation, 
jrt , A distributed parameter identification problem includes the control on the 

diffusion coefficient of the Poisson equation and a functional includes the state's 
curvature is considered. A heuristic regularization tool is introduced to manage 
codimension-one singularities during the functional analysis. Based on the 
duality principles, the approximate necessary optimality conditions is found. 
The system of optimality conditions is solved using a globalized projected 
gradient method. Numerical results, in two- and three-dimensions, implied 
the possibility of posing control on the second order derivatives and success of 
the presented numerical method. 
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wise Hessian constraint, regularization singular integral, second order control. 
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1. Introduction 



Many problems in engineering sciences and physics arc modeled by partial dif- 
ferential equations (PDEs) . The numerical solution of PDEs provides useful infor- 
mation for design engineers to predict the behaviors of corresponding systems. In 
the highly competitive world of today, it is no longer sufficient to design a system 
that performs the required task satisfactorily, but it is desired to find the optimal 
conditions. The PDE-constrained optimization (cf. [4, 14]) is an effective way to 
approach this goal. 

In the context of PDE-constrained optimal control it is common to pose either 
pointwise or global control on the state solution and/or state's derivatives. Most of 
works in this regard are focused to apply control on the the zero and/or first order 
derivatives of filed solution. 

However, most of the physical phenomena are spatial second order in nature. 
This means that the behavior of a point is a function of its neighborhood. This 
make sense to pose control on the second order derivatives of the state solution, 
in addition to the zero and/or first order derivatives. Moreover, in some cases 
we may not have sufficient information about the desired conditions in terms of 
the zero and/or first order information. Lets to give an instance, to make the 
mentioned problem more clear. In heat transfer problems, hot-centers (the locus of 
heat flux concentration) are a subset of temperature field critical points, i.e., where 
the temperature gradient approaches to zero. Now assume we want to control the 
position of hot centers. It is clear that to pose the control on the state's gradient 
does not make a physically consistent design problem, because, the cold-centers 
and saddle points are also included within the set of mentioned critical points. An 
effective way to consistently contrast between these three types of points, is to 
consider the second order information too, i.e., to include the Hessian of the state 
solution. It is worth mentioning that, some of the currently used design problems 
could be enriched/regularized by including the second order information in their 
formulation. 

However, to pose the control on the second order information is not an easy task. 
To the best of our knowledge there is neither a theoretic nor a numeric work in this 
regard. Since theoretical study on this topic is not easy using available functional 
analysis tools, we shall perform a numerical quest on this topic here. To be more 
specific, we consider a topology optimization problem (cf. [4]) in this study. Our 
problem includes the minimization of an integral functional includes the state's 
first and second order derivatives where the sate is governed by a Poisson equation. 
Moreover the diffusion coefficient of the Poisson equation is the control parame- 
ter in this study. In fact we deal with a PDE-constrained distributed parameter 
identification problem in this study. 



2. Problem statement 

Consider a simply connected domain Q C M. d (d — 2 or 3), with sufficiently regu- 
lar boundary E = dQ. where the optimization problem will take place. Consider the 
following design problem in steady-state heat (or electrical) conduction: given two 
isotropic conducting materials, with thermal conductivities a and /?, < kp < k a . 
Suppose that at each spatial point x G Q the conductivity is given in terms of the 
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a— phase characteristic function \ = x( x ): where % £ {0: 1}) as follows: 

*:(*) = Xfc« + (1 - X)kp (2.1) 

where fc Q and kp denote the thermal (or electrical) conductivity of phases a and 
(3 respectively. In fact, x(x) plays the role control or design parameter here. The 
goal is to find x( x ) configuration which solve the following design problem: 

(P) := argminF(x)= / f(x, Vu, Wu) dCl 
xeT J v 

subject to: 

V ■ (k(x)Vu) = g(x) in Q, u= u on S (2.2) 

where PCS! denotes the support domain of the integral functional, d£l denotes 
the volume measure induced in Q, u £ X U (Q) is the state function, Vu is the state 
gradient, VVm denotes the local Hessian matrix of the state function, i.e., VVu = 
[uij]dxd, i,j G {x, y} {u a b = §^jb), 9 S X g (Q) is a given source function, w (x) G 
X U (Q) denotes the Dirichlet boundary condition and T denotes the admissible 
design space which is defined as follows 

T ■■= {X e X X (Q) \ R L \Q\< f X dn<Ru\Q\, x 6 {0.1} } (2.3) 

JQ 
where < Rl < Rjj < 1 ar e given lower and upper bounds on a-phase material 
resource respectively, and |0| denotes the measure of il (volume for d = 3 and area 
for d = 2). In above formulae X u , X g and X x denotes some appropriate Banach 
spaces. Since our goal in this study is not to do a rigorous mathematical analysis, 
we simply assume that these function spaces possess sufficiently regularity required 
by our analysis. It is well known that optimal design problems like (P), are ill- 
posed and do not admit an optimal solution in class of characteristic functions (cf. 
[8]). The relaxation form of these problem is usually achieved extending the design 
domain S to the convex and continuous space 3 which is defined as follows, 

E = {weX w (Q)\ R L \Q\ < [ w rfx < Ru\Q\, < w < 1} (2.4) 

where X w is a sufficiently regular Banach space. In fact in the relaxed design 
problem the material conductivity could varies gradually between phases a and 
j3 and we have somehow a functionally graded material. The SIMP approach [4] 
is employed in this study to achieve a near 0-1 topology. In this method the 
conductivity function (2.1) is replaced by the following function: 

k(w) =w q {k a -kp) + kp (2.5) 

where q > 1 is a penalization factor which is commonly called as the SIMP power. 
In the present study, the function / in (P) is assumed to have the following form: 

/ := aw b {n- k ) c , 

where a,b,c € R, b,c > 1 are some given constants, k is the local mean curvature of 
state iso-contours, i.e., k = V ■ (Vu/|Vu|). It is clear that k is a nonlinear function 
of both the first and second order derivatives of u (cf. [15]). Therefore, this type 
of functional it is a fairly good choice for the purpose of this study. Besides it 
has a clear physical interpretation and has certain applications in practical design 
problems. For example, a negative curvature with large absolute value is signature 
of thermal centers in the heat transfer problems (cf. [18]). 
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3. REGULARIZATION OF THE CURVATURE SINGULARITIES 

Assume uo is a constant function on T, then according to the maximum princi- 
ples, there is at least a point Xo G f2 such that |Vu| =0). To avoid the curvature 
singularity in this regard, we simply replace |Vu| by |Vu| e := \Je 2 + \Vu\ 2 in the 
definition of mean curvature. Therefore the regularized mean curvature is defined 
as follows: n e — V • (Vu/|Vu| e ), where the small number e G K + is the curvature 
regularization parameter. For the purpose of convenience we simply use k instead 
of n e henceforth. In practice we use a sufficiently small value (10 -20 in this study) 
instead of e without any outer iteration on e to recover the original problem. 

Consider the discretized version of the optimal design problem with the minimum 
grid size h. According to [15], the maximum absolute value of the curvature which 
can be resolved by such a grid resolution is equal to 4 . Therefore, if the computed 
value of k locates outside of — ^ < k < i,we merely replace that value with either 
— -r or j depending on the sign of curvature. 

4. A Local analysis near a shallow gradient regions 

Suppose that u is C 2 near a generic spatial stationary point Xo, i.e., |Vu| = 
and VVu is nonsingular there. Consider d-dimcnsional closed ball B r C Q with 
radius r £ M + centered at xo . We would like to study the behavior of the objective 
functional in (P) restricted to ball B r . In this regard we can restrict ours to case 
ko = i.e., to study the behavior of the following functional: 

F\ Bt := I aw b K c dn (4.1) 

Since k is invariant under an Euclidean transformation, we assume xq = 0. Consider 
an approximation of u in the neighborhood of Xo by the Taylor expansion up to 
the second order derivatives (recall that |Vu(xo)| = 0): 

u(x) w w(x ) + - x T ■ VVu(x ) • x (4.2) 

Since n is rotationally invariant, we can assume that 

VVu(x ) = diag(A (x ), ■ ■ • , A d (xo)) (4.3) 

where A^(x ) denotes the z-th eigenvalue of VVu(x) at x = x . Therefore, by 4.2 
we have 

1 d 
u(x) Ps u(xo) + - Y. A i( x o)^ ? ( 4 - 4 ) 

8=1 

The local structure of the solution in the vicinity of xo is a function of relative values 
of Aj(xo). There are four generic types of local structures in this regard: sheet-like, 
tube-like, blob-like, double-conc-like local structures. Among these local features, 
the blob (ellipsoidal) and double-cone (hyperbolic) types have the maximum mean 
curvature concentration. For the ease of analysis, we restrict ours to the blob-like 
local structures (when Ai « • •• » A^). Assume Ai « • • • « A^ = A, i.e., a spherical 
local structure, using (4.4) and straightforward computations, we have: 
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substitution of 4.5 in 4.1 results: 

F\ B = ( aw h ( d ~ 1 = ) dfl (4.6) 

Considering the symmetric structure of the blob-like structure and using the spher- 
ical coordinate results: 

F\ Br = I aw b (— -\ n(2r) d - 1 dr = ( f r^ " 1 dr (4.7) 

where constant £ is equal to 2^ d_1 ^ (d— l) c Traw b . Considering 4.7, integral functional 
4.1 is well-defined in the vicinity of a regular shallow gradient point if c < d — 1. 
It turns out that we are not free to select the power of the curvature term in our 
integral functional. For instance, assume we are interested to solve a topology opti- 
mization problems subject to pointwise constraints on mean curvature. According 
to above analysis, in the case of two dimensions, we can not use a quadratic (or 
higher order) penalty function (cf. [5, ch. 3]) to manage local constraints (it is 
worth mentioning that using quadratic penalty functions is very common to man- 
age such problems. In the other word we have to use exact penalty functions (cf. 
[5, ch. 4]) in these cases. 

5. REGULARIZATION OF SINGULAR INTEGRALS 

To derive the first order necessary optimality conditions we have to deal with 
some singular integrals to continue the analysis. The good of this section is to 
introduce some (heuristic) regularization tools to manage such cases. Let /(x) : 
R d -> 1 be a real-valued function defined on the codimension-one C 2 manifold 
M C R d which is embedded within domain Q. We are interested to compute the 
following singular integral: 



I := / /(x) dT = / /(x) 8{M) dn (5.1) 

J M JQ 

where dT denotes the surface measure induced on .M, S(Ai) denotes the delta- 
distribution concentrated on M. Consider bjv[ G C 1 (0) as the Euclidean signed 
distance function with respect to manifold M. Assume {-)% xt denotes the smooth 
extension (along the normal direction) of scaler field (•) defined on M. into J\f(A4,a), 
where N(M.,a) denotes ±er-distance neighborhood of M.. Moreover assume that 
the manifold forms the boundaries of Af(M,a) is non-degenerate (there is no self- 
crossing). Assume that M{M. 1 cr) is spanned by sA-system such that s coincides 
to the iso-contours of the distance function &» and A coincides to normals on the 
6^4-constant surfaces. Therefore, ds is a (d — l)-dimcnsional (area) metric and 
dX is a one-dimensional metric. It is obvious that sA-system consistently spans 
N(M,cr) 1 if iso-contours of 6^1 do not cross each other (which is followed by our 
assumption). Therefore we can use sA-system together with its corresponding d- 
dimcnsional volume element, dsdX, to compute d-dimensional volume integrals. 
Applying sA-system to integral (5.1) results: 

h= f a Pa(X) f ti xt (x)dsdX (5.2) 

where lim CT ^o la = I, the manifold Mb C N{M, a) coincides to an s-snrface which 
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is located at distance b from M , and the weighting function p a E C°° (M) is chosen 
such that the intensity (measure) of /(x) be preserved, i.e., 

Pa{X)dX=l (5.3) 



\=-a 

the other requirements of function p a are: p a {x) 7^ for x E (— er, +cr), p a (x) = 
for .t ^ (— <T, +0"). In fact, p^ makes a converging sequence to the classic delta 
function, S, in the sense that: lim CT _>o Pa{x) = 5(x) (cf. [6]). in the present study, 
the following normalized Gaussian distribution is used as p a function: 

pAx) f P>- d exp(^^), if |A|<a 

where /?° is determined bases on constraint (5.3). Note that p a plays also the role 
of a regularization kernel in our formulation (cf. [7, 20]). More precisely, for every 
point xo E M., Pa does convolution of /(xo) along the line segment (— a, +cr) which 
is normal to M. and passes from xo. Therefore, if the extension of / to neighborhood 
of the desired manifold be not smooth, e.g., there are some equi-distanced points 
from the manifold surfaces, our method performs well. 

Now, lets to back from sA-coordinate to our original Cartesian coordinate. The 
coarea formula [10] is used for this purpose. Assume that the change of variable is 
global in £7, the coarea formula can be expressed as: 



|V6m| dQ= ds dX (5.5) 

'n J 

i.e., |V6.m| dil = ds dX (cf.: [9, ch. 3]). Applying (5.5) to (5.2) results: 

la = f Pa(b M (x)) /r*W \Vb M \ dn (5.6) 

JM(M,a) 

since p a = for x E Q \ N(M, a), 

Ia= [ Pa{b M {*)) /r'(x) \Vb M \ dn (5.7) 

Jn 
it is worth noticing that a formula similar to (5.7) is derived heuristically (not 
rigorously) in [15, ch. 1]. 

The Euclidean signed distance function bj^t can be efficiently computed solving 
the following eikonal equation [15, 16]: 

|V&m| = 1 in Q, b M = on M (5.8) 

this equation can be solved efficiently using either the fast marching method [16]. 
Notice that the solution of (5.8) is a member of BV(Q), i.e., it is not essentially 
sufficiently smooth. However, due to the application of the regularization kernel 
Pa in our formulation, we are not concerned about the smoothness of 6^4 . We also 
need to extend function / defined on M. along the normal direction into A/"(A / 1, a). 
The following PDE is suggested in [1] to do this job: 

f ext ■ Vb M = in Q, f ext = f on M (5.9) 

where f ext denotes the extension of / inside Af(M.,a). In [1], an efficient fast 
marching method is suggested to compute bjvi and f ext simultaneously. The im- 
plementation of this method will be used in this study. 
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For the convenience in notation, the following notation is used in this study, 
everywhere required: 



(•) dT = / (( (•) )) dfi (5.10) 

M JN{M,a) 

Let o(x) and 6(x) : R d ->• R and c(x) and d(x) : R d -> R d be functions defined 
on co-dimension one manifold A4. To do our sensitivity analysis, we need the 
following properties for operator (((•))) (note that these relations are not exactly 
hold however we are hopeful, without proof, to validity of them in the sense of 
approximation) : 

C {{ab}) dQ » f ((a)) bdQ (5.11) 

(•) •/(•) 

«c-d» dfl w / (( c )) ddn (5.12) 



6. The first order necessary optimality conditions 

There are several methods to solve a distributed parameter identification prob- 
lem. Due to large number of control parameters, a gradient based method is adapted 
in this study (cf. [3]). Lets to define the meaning of differentiability on function 
spaces. There exist several differentiability notions in mathematical programming 
literature. The notion of Gateaux derivative is applied here. 

Definition 6.1. (Gateaux derivative, cf. [3]) Consider Banach spaces Y and W 
and U as open subset of Y. A function /:[/—> W is called to be Gateaux 
diffcrcntiablc at u £ U if for every test function v £ Y the following limit exist: 

/'(«):= lim /(^)-/W 



In this case we show the Gateaux derivative symbolically by f'(u). If Y is a Hilbcrt 
space, which is assumed in this study, then f'(u) lives on the dual space of Y. 
Therefore using the Riesz representation theorem, there is a unique e £Y such that 
(e, v) = f'(u), where (•,•) denotes the inner product on Y. In this case, without 
confusion, it is common to call e as the Gateaux derivative. We use notation df(u) in 
this case, i.e., df(u) = e. It is easy to verify that under some mild conditions (which 
usually hold in practice) most of properties for classical derivatives have equivalent 
extension to Gateaux derivative. In the case of multi- variable functional, the partial 
Gateaux derivative are denoted by f',-., dr.)f(' • • ) symbols in this study. For the 
purpose of convenience the Gateaux derivative is called as the directional derivative 
in this study, henceforth. Without confusion, the notation (-, -)q := Jn(')i') dn also 
used to denote the duality pairing on function spaces in this study. 

Consider an arbitrary function p £ X P (Q), where X p is a sufficiently regular 
Banach space. Lets to introduce the following lagrangian augmenting the inner 
product of p and the Poisson equation to the objective functional in problem (P): 

C(w, u,p) ■- F(w, u) + (V ■ (kVu) - g, p) +(u- u , p) s 
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The set of points satisfy the first order necessary optimality conditions of problem 
(P), denoted by O can be expressed as follows (cf. [3]): 



O := { (w,u,p) e (H x X U (Q) x X p (Q)) 



d w C(w,u,p) = in Q (C.l) 
d u C(w,u,p) = inQ (C.2) 
d p C(w,u,p) = in Q (C.3) 



where d w C, d u C and d p C denote respectively the partial directional derivatives of C 
with respect to w, u and p along arbitrary test directions Sw G X W (Q), Su G X U (Q) 
and 5p G X P (Q) respectively. In fact, O includes constrained stationary points of 
augmented lagrangian C. Regarding to d w C in condition (C.l) of set O wc have, 

(d w C, Sw) Q = (d w f, Sw) v + (pV ■ (d w kVu), Sw) Q := h + I 2 (6.1) 

where d w f = abw t - b ^ 1 ^ (k — k ) c and d w k = qw l ^ q ~ 1 \k a — fcg). For 7 2 we have, 

h = (pd w kS7u ■ n, Sw)^ - (d w kVu ■ Vp, Sw) := I 3 + J 4 (6.2) 

where n denotes the outer unit normal on E. 

Regarding to d u C in condition (C.2) of set O wc have, 

(d u C, 5u) Q = (d u f, 5u) v + (V • (kVSu), p) Q + {Su, p> E := I 5 + 1 6 + I 7 (6.3) 

Applying two consecutive times the integration by part followed by the divergence 
theorem to term I§ in (6.3) results: 

h = {kVSu • n, p) s - (kVp ■ n, 6u) x + (V • (fcVp), <5w) Q := h + h + ha (6.4) 

Since u is fixed on S, we have Su = on S therefore /§ = Iq = in (6.4). For term 
7.5 in (6.3) we have: 

U = 'A!^, *•>„ = X 4/ V ■ (™J - 7 °7 7 " a| 3 7fa) ) «1 := A (6.5) 

To make expression concise, lets to define the following local operator P(u), 

P(u) :=l d -^L^^L (6.6) 

|Vm| |vu| 

where 1^ G R dxd denotes the identity matrix and (g> denotes the tensor product of 
ci-vectors, i.e., a®b = [aibjjdxd- Using (6.6) in (6.5), J\ can be written as follows: 

Applying the integration by part and the divergence theorem to (6.7) results: 

J BjJ VSu ■ P(u) -m^fVSu. P(u) ■ VQJ) ^ == _ 
Jap |Vu| J p |Vm| 

where m denotes the outer unit normal on dV and 91? denotes the sufficiently 
regular codimension-one manifold forms the boundaries of T>. Now lets to proceed 
the derivation, using regularization tools introduced in section 5. For J 2 in (6.8) 
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ve have 




*<*>£<< 


' <9 K / Vfc • P(«) • m\ \ JO (5.12) 
\ |Vu| // 


= [ 6u( 


7ft./P(«).mU. kdr _^ 


JdX > 


A |Vu| // .A* 


:= J4 — J5 





|Vu| 

S K /P(u) • m x 

ou ail 



\Vu\ 

(6.9) 

where X := Af(dT>, a) and k denotes the outer unit normal on boundaries of X. 
For J4 in (6.9) we have: 



jU<<Hfr=>>- k >> fc --* ,6 - 10) 



where 3^ := Af(dX, a). Applying the integration by part and the divergence theorem 
to term J 3 in (6.8) results: 

t SuV(8 K f) ■ P(u) m dr f /V(g /) ■ P(„) X u ^ == ^ _ Jg 

■/ax* |Vu| 7c V l Vu l / 

'""" u 7,(( !a ^ ))"- j - 

Collecting terms include implicit function 5u, we can solve the following adjoint 
Poisson equation to enforce the condition C.2 of set O: 

V • (fc(to)Vp) = ft(u) in Q, p = on £ (6.12) 

where w solves the direct problem (2.2), /i(u) = ft-i(w) — ^(m) — ^(u) + hi(u) and, 

ki{u) := V ' I JV«J J ^ (X) 

where Zx>(x), I^(x) and Xy(x) denote the characteristic functions of the spatial 
domains 2?, A" and y respectively, i.e., for instance 2j>(x) = 1 for all x G P and 
Xd(x) = elsewhere. According to our numerical experiments, /ii(x) <C /ii(x) 
(i = 2, 3, 4) such that hi can be ignored without a sensible loss in the accuracy of 
numerical solutions. Since p = on £, integral I3 in (6.2) is equal to zero. 

It is evident that when u solves the direct Poisson equation (2.2) the condition C.3 
holds in V. Putting altogether, the (approximate) first order necessary optimality 
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conditions for optimization problem (P) can be expressed as follows: (OC) as follows: 

V • (fc(w)Vu) = g(x) in Q 

u(x) = Wo( x ) on S 

(OC) := { V • (k(w)Vp) = h(x) in Q 

p(x) =0 on £ 

Vs(w — d w C) — w = in Q 

where d w £(x) = ji(x) - j 2 (x), ji(x) = d w f I-p(x), j 2 (x) = d w kVu ■ Vp and 
operator 'Ph(-) denotes the orthogonal projection onto the admissible control space 
S. Due to the convexity of 5 the optimality conditions (OC) is stated in terms of 
the projected gradient with respect to the admissible control space. Due to specific 
structure of S the projection 'P-(-) can be computed very efficiently in practice. 

It is worth mentioning that in the special case when T> := Q, since Su = on £ 
we have J 2 = Ji = 0. Therefor we can solve the adjoint Poisson equation exactly 
without any regularization tool and optimality system (OC) is exact in this case. 

Consider an orthogonal local curvilinear coordinate system coincides to the prin- 
cipal curvature lines on the w-constant surfaces, in fact the local (intrinsic) tangen- 
tial coordinate systems. Lets to denoted by s(x) and {ti(x)}^T 1 1 respectively the 
unit normal and tangent vector which form the basis vectors of this local curvilinear 
coordinate system. It is easy to show that (cf. [12, ch. 1]): 

d-1 

t d = s(x) + J2 *i( x ) ( 6 - 13 ) 

considering 6.6 together with 6.13, it turns out that for every d-dimensional vector 
v(x), P(u(x)) • v(x) is a <i-dimcnsional vector w(x) such that w(x) is tangential 
to u(x)-constant surfaces, i.e., P(u(x)) ■ v(x) ■ Vii(x) = 0. Therefore, at any 
Xo where V(5u(xo) be orthogonal to the u(xo)-constant surface, 8 u k, will be equal 
to zero. Therefore in a special case when dT> is defined as a it-constant surface 
hi = J2 = h 3 = and so we can solve the adjoint Poisson equation exactly without 
any regularization tool. 

7. Numerical method 

In this section we briefly mention the numerical method used to solve the opti- 
mality condition OC. The physical domain is discretized into a uniform Cartesian 
grid with an TV + 1 grid points along each spatial dimension. The direct and ad- 
joint systems are solved by the cell centered finite volume method (FVM). In this 
way, we have an N control volumes along each spatial dimension. Each control 
volume includes a degree of freedom for direct and adjoint fields ni addition to a 
degree of freedom for the control variable. All degrees of freedom are defined as 
center of control volumes. The harmonic averaging is used to compute the diffusion 
coefficient on the faces of control volumes. It is worth mentioning that, according 
to our experiments and also results reported in [11], the cell-centered FVM is free 
from the chcckcrboard-likc instability (cf. [4]) which is usually connected to finite 
element solution of topology optimization problems. 

In the present study, the spectral projected gradient method developed in [19] 
is used to solve the optimality conditions (OC). To solve the linearized optimality 
condition (OC), the value objective function together with its gradient is required 
which needs solution of direct and adjoint Poisson equations. Due to existence of 
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jump in diffusion coefficient, the condition number of the coefficient matrix related 
to the direct and adjoint problems will be very large using traditional iterative 
solvers (in particular for large values of kp/k a ). The multigrid preconditioned 
conjugate gradient method (MGCG) [17] is used to solve linear systems of equations 
in this study. The main benefits of this method arc its excellent performance in 
addition to its nearly independent convergence-rate to factor kp/k a - 

8. Results and discussion 

The success and performance of the presented approach is studied in this section 
through several numerical examples. Some of the input parameters are as follows. 
In all of the examples, the design domain is chosen to be Q = [— 0.5,0. 5] d . The 
following strict equality resource constraint is considered Rl = Ru — 0.5 and 
the threshold for the projection onto the admissible control domain is taken equal 
to l.e — 6. To keep the initial design inside the feasible domain, we start with 
w(x) = 0.5. The smoothing width of the regularization kernel, a, is equal to 3Aa;, 
where Ax is the grid spacing. The convergence threshold of MGCG algorithm is 
l.e — 20, except where it is given explicitly. The jacobi method is used as our 
multigrid smoother and the number of smoothing iteration on each grid level (after 
and before the prolongation and restriction operations) is one. We consider two 
form of definitions for the support of integral in the objective functional (T>). In 
the first form, T> is defined by a spatial function and in the second form, T> is taken 
to be a function of the field solution, for instance, we first solve the direct problem 
and then patch a portion of the spatial domain based on the computed local mean 
curvature value. In fact, this type of definition can be considered as an ingredient 
of our ultimate goal which is the pointwise control on the second order derivatives. 
In this context, the present works plays the role of the sub-problem solver for an 
augmented lagrangian method (cf. [13]). The optimization cycle is stopped when 
the difference between two consecutive topologies be below 0.1%. All floating point 
arithmetic is performed to the double precision accuracy. A personal computer 
with an AMD 2.41 GHz CPU and 2.5GB RAM is used as the computing platform. 

8.1. Two dimensional results. The following examples are used to evaluate the 
presented method in two spatial dimensions. In these examples, the physical domain 
is divided into a 256 x 256 uniform Cartesian grid. 

Example 8.1. k a = 2, hp = 1, q = 1, <?(x) = 1, uq(x) = 0, a = —1, b = 0, c = 1, 

k q = 0, V = { x 6 [-0.5, 0.5] 2 | |x| < 0.25, \y\ < 0.25}, where x = (x, y) T . 

Example 8.2. this example is like to 8.1, else 6=1. 

Example 8.3. this example is like to 8.1, else the conductivity ratio is increased 
to 200, i.e., k a — 200, kp = 1, the SIMP power is also increased to q = 5 avoid 
intermediate densities. 

Example 8.4. this example is like to 8.1, else the objective function domain is 
defined as a function of the state curvature at the start of the optimization, i.e., 
assuming kq = —6 the characteristic function x('Z') is defined as: 

1 if k(x) < Kq 



x(2?) ' if k(x) > k 
Example 8.5. this example is like to 8.4, else 6=1. 
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Example 8.6. this example is like to 8.4, else the direction of the optimization is 
reversed, i.e., a = 1. 

Example 8.7. this example is like to 8.4, else the conductivity ratio is increased 
to 200, i.e., k a = 200, kp = 1, and the SIMP power is increased to 5 too. 

Example 8.8. this example is like to 8.7, else the direction of the optimization is 
reversed, i.e., a = 1. 

Example 8.9. this example is like to 8.4, else the right hand side of the direct 
problem is changed to: g(x) = cos(37ra;) cos(37ry). 

Example 8.10. this example is like to 8.9, else 6=1. 

Example 8.11. this example is like to 8.9, else the direction of the optimization 
is reversed, i.e., a = 1; also q = 10. 

Example 8.12. this example is like to 8.11, else 6=1. 

Example 8.13. this example is like to 8.4, else the asymmetric right hand side 
<?(x) = cos(7ra:) sin(7rj/) is applied. 

Example 8.14. this example is like to 8.13, else 6=1. 

The resulted topologies related to examples 8.1- 8.14 are shown in Figures 1- 
3. Figures 4-6 show the variation of the objective functional as the optimization 
proceeds. Plots clearly show that, the presented method is enable to effectively 
reduce the objective functional in all cases and to move toward the optimal solution 
(if there be any). Moreover, in most of cases a near 0-1 topology is achieved. 
Therefore, we can conjecture about the existence of an optimal solution for problem 
(P) defined in this study. As it is clear from the plots, we do not have essentially a 
monotonic reduction in the objective functional value which is an expected result 
due to the use of a nonmonotonic line-search in our optimization algorithm (cf. 
[19]). 

To study the stability and convergence of our method with respect to the grid 
refinement, example 8.4 is considered with grid resolutions: 2™ x 2", n = 5, ■ • • 10. 
The resulted topologies related to this numerical experiment is shown in Figure 
7. Plot shows that the macroscopic features of the topology have an excellent 
convergence. Moreover, the microscopic features follows the same trend without 
topological instability. Therefore, we can conjecture on the well-posedness and 
stability of the applied numerical method. Notice that, a micro-scale topological 
convergence is generally not expected due to this well-known fact that by the grid 
refinement the optimal solutions tend to form smaller and smaller microstructures 
(cf. [2]). 

To make sense about the direct and adjoint fields at optimal solutions, the direct 
and adjoint fields at the optimal solution corresponding to examples 8.4 and 8.8 
are plotted in Figure 8. 
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Figure 1. Optimal design in two spatial dimensions: a-f are re- 
lated to the final topologies of examples 8.1-8.6 respectively. 
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Figure 2. Optimal design in two spatial dimensions: a-f are re- 
lated to the final topologies of examples 8.7-8.12 respectively. 
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Figure 3. Optimal design in two spatial dimensions: a and b are 
related to the final topologies of examples 8.13 and 8.14 respec- 
tively. 



Figure 4. Optimal design in two spatial dimensions: the objective 
function history, a-f are related to examples 8.1-8.6 respectively. 
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Figure 5. Optimal design in two spatial dimensions: the objective 
function history, a-f are related to examples 8.7-8.12 respectively. 
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Figure 6. Optimal design in two spatial dimensions: the objective 
function history, a and b are related to examples 8.13 and 8.14 
respectively 
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Figure 7. Optimal design in two spatial dimensions: the grid 
resolution study, a-f are related to the final topologies computed 
on the grid resolutions 2" X 2™, n = 5, ■ • • , 10 respectively. 
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Figure 8. Optimal design in two spatial dimensions: the direct 
(left) and adjoint (right) fields at the final design for examples 8.4 
(top) and 8.8 bottom. 



To make sense about the computational cost of the presented approach, example 
8.4 is considered with different conductivity ratios jj*- = 2, f 0, 100, and grid reso- 
lutions: 2" x 2", n = 5, • • • 11. Notice that for the conductivity ratio 100, SIMP 
power 5 is considered. 

Results of this numerical experiment are shown in Table 1. Table illustrates that 
the number of optimization cycles has a little dependency on the grid resolution and 
the conductivity ratio. The number of MGCG iterations is increased by increasing 
the conductivity ratio, r^. But for a fixed r^, the number of MGCG iterations 
is almost independent from the grid resolution, which show the reliability of this 
solver for large-scale problems. The reported computational cost illustrates the 
excellent performance of the applied numerical strategy in this study. To the best 
of our knowledge such a large number of design parameters (2 22 ) is not reported in 
literature of topology optimization so far. 



8.2. Three dimensional results. In this section we show the success of the pre- 
sented method in three spatial dimensions. For this purpose the following example 
is considered on grid resolutions: 2™ x 2" x 2", n = 5, 6, 7. 
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Table 1. Performance analysis: results include total optimization 
iterations (itero), total functional evaluation (nfunc), total gradi- 
ent evaluation (ngard), total MGCG iterations (iterm) and total 
computational cost in second (cpu); for different conductivity ra- 
tios (r/j = k a /kp) and different grid resolutions (grid). 



grid 


rk 


itero 


nfunc 


ngrad 


iterm 


cpu (sec) 


(2 6 ) 2 


2 


10 


11 


11 


462 


0.25 


(2 7 ) 2 


2 


10 


11 


11 


484 


2.15 


(2 8 ) 2 


2 


10 


11 


11 


486 


10.62 


(2 9 ) 2 


2 


11 


12 


12 


528 


46.76 


(2 10 ) 2 


2 


11 


12 


12 


542 


203.1 


(2 11 ) 2 


2 


13 


12 


12 


612 


1180.6 


(2 6 ) 2 


10 


10 


11 


11 


792 


0.36 


(2 7 ) 2 


10 


10 


11 


11 


748 


3.16 


(2 8 ) 2 


10 


10 


11 


11 


836 


16.25 


(2^)2 


10 


11 


12 


12 


990 


78.30 


(2 10 ) 2 


10 


11 


12 


12 


1020 


382.4 


(2 11 ) 2 


10 


14 


15 


15 


1172 


1905.1 


(2 6 ) 2 


100 


10 


11 


11 


1430 


0.56 


(2 7 )2 


100 


10 


11 


11 


1716 


6.49 


(2 8 ) 2 


100 


10 


11 


11 


1958 


35.05 


(2») 2 


100 


12 


13 


13 


3014 


220.7 


(2 10 ) 2 


100 


12 


13 


13 


4940 


1573.5 


(2 11 ) 2 


100 


14 


15 


15 


6582 


9941.1 



Example 8.15. k a = 2, kp = 1, q = 1, g(x) = 1, uo(x) = 0, a = —1, b = 0, c = 1, 
and assuming ko = — 3, the characteristic function x(^*) is defined as 

1 if k £ (x) < K 



X ^°> \ if k £ (x) > K 

Figures 9 and 10 show the final topologies and objective function history cor- 
responding to three dimensional numerical examples. In all cases the number of 
optimization iterations was below 13 and the number of objective functional and 
its gradient evaluation was equal. Plots illustrate the convergence and stability of 
the presented method in three spatial dimensions too. 
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Figure 9. Optimal design in three spatial dimensions: a-c arc 
related to the final topologies computed on the grid resolutions 
2™ x 2™ x 2™, n = 5, 6, 7 respectively. 
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Figure 10. Optimal design in three spatial dimensions: the ob- 
jective function history, a and b are related to grid resolutions 
2 6 x 2 6 x 2 6 and 2 7 x 2 7 x 2 7 respectively. 



9. Summary 

The idea of posing control on the second order derivatives of the PDEs solutions 
is introduced in this study. To be specific, a topology optimization problem corre- 
sponding to the Poisson equation in which the objective functional is a nonlinear 
function of first and second order field derivatives is considered. Introducing some 
functional regularization tools, the (approximate) first order necessary optimality 
conditions is derived and is numerical solved using an appropriate gradient descent 
method. Numerical results are provided for a variety of test cases in two and three 
spatial dimensions. The Numerical results in this study, confirms the stability, con- 
vergence and efficiency of the presented approach. In fact, the numerical results 
answer to our quest in this study; the possibility of posing some degrees of control 
on the second order derivatives. 
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